#ifndef _FUNCTION_H_
#define _FUNCTION_H_

#include <iostream>
#include <limits>
#include <cmath>
#include <algorithm>
#include <vector>
using namespace std;

class Function
{
public:
    virtual double operator()(double _x) = 0;
    virtual double diff1(double _x) {
        return 0;
    }
    virtual double diff2(double _x){
        return 0;
    }
};

class Polynomial : public Function
{
private:
    double x0;
    int n;
    vector<double> A;
public:
    Polynomial(vector<double> _A, double _x0): A(_A), x0(_x0) {}
    double operator()(double _x){
        double y = A[0];
        n = A.size();
        for (int i = 1; i < n; i++){
            double dy = A[i];
            for (int j = 0; j < i; j++)
                dy = dy * (_x - x0);
            y = y + dy;
        }
        return y;
    }
};

class Discrete : public Function
{
private:
    vector<double> X, Y, dX1, dX2;
    int n;
public:
    Discrete(vector<double> _X, vector<double> _Y, vector<double> _dX1, vector<double> _dX2): X(_X), Y(_Y), dX1(_dX1), dX2(_dX2) 
    {
        if ((X.size() != Y.size()) || (Y.size() != dX1.size()) || (dX1.size() != dX2.size())){
            cout<< "Error!" <<endl;
            exit(0);
        }
        n = X.size();
    }
    double operator()(double _x)
    {
        for (int i = 0; i < n; i++){
            if (X[i] == _x)
                return Y[i]; 
        }
        return 0;
    }
    double diff1(double _x)
    {
        for (int i = 0; i < n; i++){
            if (X[i] == _x)
                return dX1[i]; 
        }
        return 0;
    }
    double diff2(double _x)
    {
        for (int i = 0; i < n; i++){
            if (X[i] == _x)
                return dX2[i]; 
        }
        return 0;
    }
};

class Bsplines : public Function
{
private:
    int dimention;
    vector<double> t;
public:
    Bsplines(int _dimention, vector<double> _t): dimention(_dimention), t(_t)
    {
        if (dimention != t.size() - 2)
	    {
		    cout<< "Error!" <<endl;
            exit(0);
        }
    } 
    double operator()(double _x)
    {
        vector<double> t1(t.begin(), t.end() - 1);
        vector<double> t2(t.begin() + 1, t.end());
        Bsplines B1(dimention - 1, t1);
        Bsplines B2(dimention - 1, t2);
        //Definition 3.21
        if (dimention == 1)
        {
            if (_x > t[0] && _x < t[1]){
                return (_x - t[0]) / (t[1] - t[0]);
            }
            else if (_x >= t[1] && _x < t[2]){
                return (t[2] - _x) / (t[2] - t[1]);
            }
            else{
                return 0;
            }
        }
        else if(dimention > 1)
        {
            return (_x - t[0]) * B1(_x) / (t[dimention] - t[0]) + (t[dimention + 1] - _x) * B2(_x) / (t[dimention + 1] - t[1]);
        }
    }
};

#endif
